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Using the Constrained Path Monte Carlo (CPMC) 
method, we simulated the two-dimensional, three-band Hub- 
bard model to study pairing, charge, and spin correlations 
as a function of electron and hole doping and the Coulomb 
repulsion Vpd between charges on neighboring Cu and O lat- 
tice sites. As a function of distance, both the d^2_y2-wave 
and extended s-wave pairing correlations decayed quickly. In 
the charge-transfer regime, increasing Vpd decreased the long- 
range part of the correlation functions in both channels, while 
in the mixed- valent regime, it increased the long-range part 
of the s-wave behavior but decreased that of the d-wave be- 
havior. Still the d-wave behavior dominated. At a given 
doping, increasing Vpd increased the spin-spin correlations in 
the charge-transfer regime but decreased them in the mixed- 
valent regime. Also increasing Vpd suppressed the charge- 
charge correlations between neighboring Cu and O sites. Elec- 
tron and hole doping away from half-filling was accompanied 
by a rapid suppression of anti- ferromagnetic correlations. 



I. INTRODUCTION 

In this paper, we will report the results of a quantum 
Monte Carlo (QMC) study of the ground-state proper- 
ties of the two-dimensional, three-band Hubbard model. 
Of the three simple electronic models commonly stud- 
ied as possible models of the cuprate superconducting 
materials, the three-band Hubbard model has been the 
least intensively studied, partially because of the general 
belief that its low energy excitation spectrum is similar 
to the other two. Indeed in the strong coupling limit, 
both the one-band and three-band Hubbard models have 
the t-J model as an approximate limit. However, there 
still remains some controversy about whether one-band 
models, like the Hubbard and t-J models, are adequate to 
describe the low energy physical properties of the cuprate 
superconductors. One of our objectives was to study the 
possible existence of superconductivity in the three-band 
model in regions where model parameters are physical as 
opposed to regions where asymptotic models are clearly 
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more appropriate. A focus of our ground state study is 
the effect of the inclusion of Vpd, a repulsive Coulomb 
interaction between charges on neighboring Cu and O 
lattice sites. 

The most solid information about possible supercon- 
ductivity in the one-band Hubbard model has come from 
a series of QMC calculations. For instance, riising a fi- 
nite temperature QMC method. White et al.El found an 
attractive effective pairing interaction in the d,j,-2_^ and 
extended s-wave channels. Moreo and Scalapinocl sub- 
sequently found pairing correlations but also found that 
they did not increase when the lattice size was increased. 
This result, suggesting the absence of off-diagonal long 
range order, was consistent with an earlier QMC study 
by Imada and Hatsugaifl 

The fermion sign problem limits almost all QMC cal- 
culations of the Hubbard model to small system sizes 
and finite temperature simulations to high temperatures. 
These limitations had left open the possibility that super- 
conductivity still lurked at larger systems sizes and lower 
temperatures. Recently, a new zero temperature QMC 
method,|,the Constrained Path Monte Carlo (CPMC) 
method,Q was developed to get jaround the sign problem. 
Using this method, Zhang et aljj calculated ground-state 
pairing correlation functions for the Hubbard model as a 
function of the distance and found that as the system size 
or the interaction strength is increased the magnitude of 
the long-range part of the pairing correlation functions 
vanished for both the dx^_y2 and extended s-wave chan- 
nels. Although this method produces an approximate 
solution, its results, together with the null results from 
previous QMC studies, are very discouraging for finding 
superconductivity in the one-band Hubbard model. 

At finite temperature, past numerical work on the 
three-band ijwlel has included several sets of QMC 
simulations J3~0 focusing on the magnetic superconduct- 
ing and insulating properties of the model. As for the 
QMC simulations of the Hubbard model, the sign prob- 
lem limited these studies to relatively high temperatures 
and small systems. In fact the sign problem for the three- 
band model is probably more severe than for the one 
band model. At least the one band model lacks a sign 
problem at half filling. In general, an anti-ferromagnetic 
state is found at half-filling which is strongly suppressed 
upon doping. Attractive interactions between pairs were 
found with a spectrum of results: Dominance of ex- 
tended s-wave and d-wave pairing have been separately 



reported, leading to claims of extended s-wave and d- 
wave ODLRO. These results remain controversial. 

At zero temperature, past numerical work on the 
three-band model has maipljy j-eonsisted of exact di- 
agonalization computations,t2l^EJ where calculations of 
hole binding eBjetgi^ was emphasized, and several QMC 
computationsE3ll30 where pairing calculations were em- 
phasized. The exact diagonalization studies unequivo- 
cally established that holes can bind. The QMC stud- 
ies, established the existence of an extended s-wave and 
d^2_y2 wave attractive pairing interaction, with one claim 
of of no evidence of s-wave superconductivity. More rjfe. 
cently, a CPMC study by Guerrero and Gubernatist^ 
confirmed the exact diagonalization result that holes bind 
but found that increasing system size tended to decrease 
the long range part of the pairing correlations. In con- 
trast to the exact diagonalization work, this QMC study 
found hole binding in the absence of a Coulomb repulsion 
Vpd between charge on neighboring Cu and O sites. The 
exact diagonalization studies found that hole binding re- 
quired an unphysically large value of Vpd- 

In the work reported here, we applied the CPMC 
method to the three-band Hubbard model and computed 
the static pairing, charge, and spin correlation functions 
for systems with 6x6 unit cells. We set the hopping and 
on-site Coulomb parameters as expected physically rele- 
vant values and studied the properties of the systems as a 
function of charge-transfer energy, electron and hole dop- 
ing, and Vpd through a range of physically relevant values. 
The consequences of varying these parameters for the 
most part depended on whether the value of the charge- 
transfer energy placed the model in a charge-transfer or 
mixed valent regime. The cuprate superconductors are 
believed to be in the charge-transfer regime. 

The remainder of our report is organized as follows: In 
Section y, we define the Hamiltonian and the physical 
quantities calculated and discuss the choice of model pa- 
rameters. In Section [II we briefly describe the CPMC 
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method, and then in Section IV, we present our numeri- 
cal results. Finally in Section^, we discuss in detail our 
main conclusions. 



II. THREE-BAND HUBBARD HAMILTONIAN 

As proposed by Emeryjij the three-band model mimics 
the Cu02 layer in the cuprate superconductors by hav- 
ing one Cu and two O atoms per unit cell, with the Cu 
atoms arranged on a square lattice and the O atoms cen- 
tered on the edges of the square unit cells. In this layer, 
Emery assumed that the relevant orbitals are just those 
of copper 3da:2„j,2 and oxygen 2px and 2py. 

The Hamiltonian has the form: 

^= E *^p(p>..+pU-)+^pE"'^ + t^pE4 
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In writing the Hamiltonian, we adopted the convention 
that the operator d]^ creates a hole with spin cr at a 
Cu 3(ij,2_j,2 orbital and p!;^ creates a hole with spin a in 
an O 2px or 2py orbital. Ud and Up are the Coulomb 
repulsions at the Cu and O sites, €d and Cp are the corre- 
sponding orbital energies, and Vpd is the nearest neighbor 
Coulomb repulsion. As written, the model has a Cu-0 
hybridization t'j'^ = Upd with the minus sign occurring 
for j = i + x/2 and j = i — y/2 and also hybridizaton 
tjjl = ±tpp between oxygen sites with the minus sign oc- 
curring for k = j — x/2 — y/2 and k = j + x/2 + y/2. 
These phase conventions are illustrated in Fig. 1. 

The values of the parameters in the Hamiltonian have 
been estimated by a number of different constrainndjien- 
sity functional and quantum cluster calculations.Ej~tJ In 
electron-volts reasonable ranges for these values seem to 
be: tpd = 1.3 - 1.6, Ud = 8.5 - 10.5, e = ep - Cd = 3.6, 
Up == 4.0-7.5, tpp = 0.65, and Vpd = 0.6-1.2. Taken to- 
gether, these estimates dcflne a reasonably limited range 
of the parameters for the which the model might be la- 
beled as "physical." 

The Cu site Coulomb repulsion Ud is a large energy, 
making doubly occupancy of Cu sites by two holes very 
unfavorable. The next largest parameter, the charge 
transfer energy e = Cp — e^ > 0, plays a special role. De- 
pendent on the relative values of e, Ud and the bandwidth. 
W, the system can be classifled in different regimes:t3 
the charge-transfer regime with Ud > e > W or the 
mixed- valent regime with Ud > W > e, where W is 
some measure of the width of the lower band. Estimates 
place the cuprate superconductors in the charge-transfer 
regime. The role of the Cu-0 hybridization tpd is im- 
portant. Through the super-exchange mechanism, this 
hybdization generates an antifcrromagnetic exchange in- 
teraction between the spins in the Cu sites. The 0-0 
hybridization tpp and the O site Coulomb repulsion Up 
are the two smallest energies. When tpp is non-zero, the 
non-interacting band structure has features that seem to 
appear in the normal state properties of the cuprate ma- 
terials. We were mainly interested is studying the conse- 
quences of Vpd- These consequences were studied in part 
by previous QMC simulations of Dopf et alB and Scalet- 
tar et alE In what follows we will scale all the energies 
by tpd- 

For the non- interacting case {Up = Up = Vpd = 0), the 
band structure is easily determined numerically and is 
illustrated in Fig. 2 for a set of parameters used. We will 
add and remove holes from lower band, and this band 
is said to be half-filled when there is one hole per unit 
cell. At half-filling, for a wide range of parameters, the 
ground-state is an anti-ferromagnetic insulator just like 
^"■jlthe one-band Hubbard model. This band gap can be esti- 
mated from Fig. 2. A better physical feel can be obtained 
/i\from the exact expression easily obtained if tpp = 0: It 
is simply e. The band-width W = ^(e/2)2 + 8 - e/2. It 



varies monotonically from 8 to zero as e varies from zero 
to infinity, e = 2 marks a value for which W = e. This 
picture does not change much for relatively|-,small non- 
zero values of tpp. As shown by Dopf et al.Ll, for e = 1 
the charge transfer gap is vanishing small, whereas for 
e = 3 a finite charge transfer gap arises in the strong- 
coupling region, e = 1 is in the mixed valent regime, 
while e = 3 is in the charge transfer regime. And we 
will see that e = 2 behaves more like the charge-transfer 
than the mixed-valent regime. In Fig. 3, we show the 
Fermi surfaces for infinite systems at the various dopings 
studied. 

li e '^ Ud, the three-band model maps into a one-band 
model with tejf ~ ^prf/f and U = Ud- For Ud ^ teff, 
the one-band model can in turn be mapped pto the t-J 
model with J = At'^fs/Ud- Zhang and Rica23 have ar- 
gued that the t-J model can also be appropriate when 
< tpd ^ e, Ud, Ud — e. Ja real materials, e/tpd is es- 
timated to be ~ 2.7 — 3.7.Efl Therefore, besides the lack 
of conclusive evidence that the one-band model super- 
conducts, it is also unclear that the mapping among the 
most studied models is appropriate for physical values of 
the parameters. 

With the numerical method used, a variety of expecta- 
tion values can be computed. We focused on the pairing, 
spin, and charge correlation functions. More specifically 
we computed the extended s-wave and the d^2_y2 pairing 
correlations as functions of distance 
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with 5 — ix, ±y. For the extended s-wave pairing 
fs*{5) = 1 for all 5 and for the d^^-y^ pairing, fd{S) = 1 
for S = ±x and fd{S) = — 1 for (5 = ±y. The magnitude of 
these quantities are dominated by a large peak in P^ (R) 
when i? = |_R| is less than a few nearest neighbor dis- 
tances. Over these distances. Pa measures local correla- 
tions among spin and charge, has little information about 
long-range pairing correlations, and may give a "false 
positive" indication of enhanced pairing. Because of this 
we will report neither the q — spatial Fourier transfor- 
mation nor the partial sums like SrJL) = X]fl<L -^aiR) 
as done in some previous works Jj^Ej Instead we will re- 
port Sa{L) — J2r>l Pa{R) where L is about two lattice 
spacings. We will also report the "vertex contribution" 
tO|ihe correlation functions (see, for example. White ct. 
al.ai defined as follows: 



where Pa{R) is the contribution of dressed non- 
interacting propagator: for each term in Pa{R) of the 
form < cl Cfc!c| >, Pa{R) has a term like < clc| >< 

clc| >. We found that in most cases the conclusions 
remain the same no matter which quantity we look at. 

For the static spin-spin correlation function we used 
the Fourier transform of the spin-spin correlation func- 
tion for the spin on the Cu sites 



(4) 



where / and m refers to the Cu sites and N is the num- 
ber of unit cells. For charge-charge correlations we com- 
puted a charge-transfer correlation function involving the 
O sites neighboring a Cu site quantity: 



C{k) 



l-Y^e^<^'^Hp{i)pm. 



(5) 



where j are the nearest-neighbors of i, p{i) = nf — n^" — 
n^^ with nf, nf^, and n^^ being the charge-density oper- 
ators on the Cu, x-axis O, and y-axis O in the unit cell 
i. 



III. NUMERICAL METHOD 

Our numerical method, the constrained path Monte 
Carlo (CPMC) method,, is extensively described and 
benchmarked elsewhere&cl. Here we only discuss its ba- 
sic strategy and approximation. In the CPMC method, 
the ground-state wave function \^q) is projected from a 
known initial wave function |V'r) by a branching random 
walk in an over-complete space of Slater determinants 
\4>). In such a space, we can write |V'o) = '^d>x{4>)\4>)- 
The random walk produces an ensemble of |0), called ran- 
dom walkers, which represent \'iPq) in the sense that their 
distribution is a Monte Carlo sampling of x{^): that is, a 
sampling of the ground-state wave function. More specif- 
ically, starting with some trial state |'0t), we project out 
the ground state by iterating 
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where Et is some guess of the ground-state energy. Pur- 
posely At is a small parameter so for if = T -|- 1/ we can 
write 



g-Arff _ g-ArT/2g-Aryg-ArT/2 



(7) 



where T and V are the kinetic and potential energies. 

For the study at hand, the initial state |V^t) is the 
direct product of two spin Slater determinants, i.e.. 



Va{R) ^ PaiR) ~ Pc.{R) 



(3) 



iV'T)=ni 



(8) 



Because the kinetic energy is a quadratic form in the cre- 
ation and destruction operators for each spin, the action 
of its exponential on the trial state is simply to transform 
one direct product of Slater determinants into another. 
While the potential energy is not a quadratic form in 
the creation and destruction operators, its exponential 
is replaced by sum of exponentials of such forms via the 
discrete Hubbard-Stratonovich transformation. For the 
on-site Coulomb term, this transformation is 
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nearest neighbor Coulomb repulsion term, we make the 
same type of transformation but we have to do it many 



more times: nfn 
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nf.nP^ + nf^nl^ + nf^n^^ + nf^n^^. 



For each term and each j in the x and ij directions, 
a Hubbard-Stratonovich transformation is required for 
a total of 8 such transformations. Because the com- 
putational time scales with the number of Hubbard- 
Stratonovich transformations, having a Vpd 7^ increases 
the computational cost by a factor of 8. 

One consequence of the Hubbard-Stratonvich transfor- 
mation is the factorization of the projection into an up 
and down spin part. Accordingly we re-express the iter- 
ation step as 



ll\cj,'J^JdxP{x)l[B, 
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(10) 



where x — {xi,X2, ■ ■ ■ ,xn) is the set of Hubbard- 
Stratonovich fields (one for each lattice site), N is the 
number of lattice sites, P{x) = (5)^ is the probability 
distribution for these fields, and B„{x) is an operator 
function of these fields formed from the product of the 
exponentials of the kinetic and potential energies. 

The Monte Carlo method is used to perform the multi- 
dimensional integration over the Hubbard-Stratonovich 
fields. It does so by generating a set of random walkers 
initialized by replicating \^t) many times. Each walker 
is then propagated independently by sampling a x from 
P{x) and propagating it with B{x). After the propaga- 
tion has "equilibrated," the sum over the walkers pro- 
vides an estimate of the ground-state wave function \iPq). 

We used two different estimators for the expectation 
values of some observable O. One is the mixed estimator 
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and the other is the back-propagated estimator 
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where ji/io) is the QMC estimate of the ground state and 
£ is typically in the range of 20 to 40. For observables 



that commute with the Hamiltonian, the mixed estima- 
tor is a very accurate one and converges to the exact 
answer as |V'o) converges to exact ground state. For ob- 
servables that do not commute with the Hamiltonian, 
like correlation functions, the back-propagated estima- 
tor has been found to give very accurate estimates of 
ground-state properties. Significant differences between 
the predictions of these two estimators often exist. 

To completely specify the ground-state wave function 
for a system of interacting electrons, only determinants 
satisfying (V'ol'/') > are needed because |V'o) resides in 
either of two degenerate halves of the Slater determinant 
space, separated by a nodal surface N that is defined by 
(■0010) = 0. The degeneracy is a consequence of both 
|0o) a-iid — IV'o) satisfying Schrodinger's equation. The 
sign problem occurs because walkers can cross N as their 
orbitals evolve continuously in the random walk. Asymp- 
totically they populate the two halves equally, leading to 
an ensemble that has zero overlap with |'0o)- If N were 
known, we would simply constrain the random walk to 
one half of the space and obtain an exact solution of 
Schrodinger's equation. In the constrained-path QMC 
method, without a priori knowledge of N, we use a trial 
wave function |V't) and require (V'tI'/') > 0. This is what 
is called the constrained-path approximation. 

All the calculations reported here were done for 
copper-oxide planes with periodic boundary conditions. 
Mostly, we study closed shell cases, for which the cor- 
responding free-electron wave function is non-degenerate 
and translationally invariant. For 6x6 unit cells, the 
dopings, producing Fermi surfaces in Fig. 2, correspond 
to closed shell fillings. In these cases, the free-electron 
wave function, represented by a single Slater determi- 
nant, is used as the trial wave function \iIit)- The use 
of an unrestricted Hartree-Fock wave function as |V't) 
generally produced no significant improvement in the re- 
sults. At half-filing, which is not a closed shell case, we 
used a linear combination of two degenerate Q = (tt, tt) 
spin-density wave states. 

In a typical run, the average number of walkers was 
600, and the time step was 0.03. We performed 1600 
steps before we started taking measurements, and we 
did the measurements in 30 blocks of 320 steps each to 
ensure statistical independence. Back propagation mea- 
surements had 40 backward steps. 



IV. RESULTS 

As mentioned before, all our simulations were done 
on lattices of 6 x 6 unit cells. For this size, 36 holes 
corresponds to a half-filed case. In units of tpd, we set 
Ud — Q, Up — Q, and tpp ~ 0.3 for most studied cases. 
We varied Vpd between and 1 for several different hole 
fillings and values of the charge-transfer energy e. We 
were mainly concerned with hole doped cases, number of 
holes. 



A. Charge Correlation Functions 



B. Spin Correlation Functions 



In Fig. 4 we show the expectation values of the charge 
on the Cu sites as a function of Vpd for several band fill- 
ings and values of the charge-transfer energy e. When 
Vpd = 0, we see that for the half- filled case, even with 
a relatively large value of e, there are substantial holes 
distributed on the O sites. When doped to 42 and 46 
holes, most of the added holes go to the O site. Except 
for the e = 1 case, increasing Vpd from transfers some 
of the O charge to the Cu. The transfer rate increases if 
e is increased. These are expected results: in the charge- 
transfer regime, Ud > e > W, with a repulsive Vpd, it 
become energetically favorable for some charge to move 
from O to Cu even at the expense of some unfavorable 
double occupancy of the Cu site caused by a large Ud- On 
the other hand, in the mixed-valent regime, Ud > W > e, 
when e = 1, we see a movement of holes from the Cu 
site to the O sites. Here the energy difference between 
the Cu and O states is smaller, and a strong on-site re- 
pulsive Ud favoring charge removable from the Cu sites 
dominates the smaller repulsive Vpd, opposing the move- 
ment of charge to the O sites. In general, the presence of 
Vpd seems to expand the charge-transfer regime. Similar 
results have been seen in finite-temper|ature QMCEI and 
zero-temperature exact diagonalizatioru studies. 

Another effect of increasing Vpd is the decreasing of the 
correlation between the charge on the Cu and the neigh- 
boring O sites. This is shown in Fig. 5. At a given band 
filling, increasing e only has a little effect. Similar be- 
havior was also seen in a finite-temperature QMC study. 
One also observes that in the charge-transfer regime the 
decreasing rate seems independent of the filling and e. 
Because C(0) « Nh/N — (ncuno), it is no surprise to 
observe an increase in C(k = 0) with an increasing Vpd- 

It is pstructive to compare our findings with Stephan 
et al.'aij results. Their exact diagonalization results 
showed that when Ucu = Uq = oo, doped holes make the 
hole on the Cu sites transfer to the O sites, and with in- 
creasing Vpd, the charge on the Cu site decreases continu- 
ously. They also found that the addition of a second hole 
produces a smaller additional change in the neighboring 
Cu-0 charge correlation than the first one does. Both be- 
haviors of the hole on the Cu site and neighboring Cu-0 
charge correlation indicate that a charge-transfer "bipo- 
laron" forms in this system. Therefore they concluded 
that the binding energy is obtained from electronic po- 
larization. From our simulation results, in the charge 
transfer region, the charge on the Cu site increases with 
hole doping and increasing Vpd, suggesting that electronic 
polarization has no effect in the physically relevant re- 
gion. Hence we expect that the binding energy is mainly 
gained from magnetic mechanism. 



In Fig. 6 we show the behavior of the local magnetic 
moment on the Cu sites. First, we see that increasing the 
hole doping increases the moment. More specifically, at a 
given doping, increasing e increases the moment. When 
e > 1, increasing Vpd increases the moment, but when 
e — 1, increasing Vpd decreases it. Clearly, the increase 
of the moment is strongly correlated with the increase 
and decrease of charges on the Cu sites. 

In Figs. 7 and 8 we show the wavevector dependence of 
the Fourier transform of the static spin-spin correlation 
function for the Cu sites as a function of Vpd for a doping 
to 42 holes. This function is plotted along high symmetry 
lines in the first Brillouin zone. The different figures cor- 
responds to different values of e. In mixed valent regime, 
Fig. 7, we see that increasing Vpd suppress this function 
over the entire zone. This suppression is consistent with 
the suppression of the local moment seen in Fig. 6. On 
the other hand, in the charge-transfer regime, Fig. 8, in- 
creasing Vpd enhances this function, again consistent with 
the enhancement of the local moment seen in Fig. 6. 
By comparing the two figures we see that for a given 
Vpd increasing e increases these correlations and sharp- 
ens the peaks in the functions. In each figure there are 
two principal peaks: One is connected with the displace- 
ment of the antiferromagnetic peak to (tTjTt — i5). The 
other is the appearance of an incommensurate structure 
at (tt ~ 5' ,n — S'). A weaker spin-density wave structure 
is at (7r,0). 

Prepintts QMC simulations of the one-band Hubbard 
modelEa'EZl have seen a similar shifting (and splitting) of 
the peak in the static structure factor from the antifer- 
romagnetic position (tt, tt) to positions (tt, tt — d) on the 
face of the Brillouin zone and (tt — S',7t — S') along the 
diagonal direction. This behavior is in agreement with 
the experimental data for LSCO presented in Fig. 3 of 
Ref.Ea, where a minimum is observed at (7r,7r) along the 
diagonal direction. Our CPMC simulations of < — t — U 
Hubbard model also found that for a large t , a weak 
peak appears along the diagonal direction. We remark 
that we did not study the dependence of either of these 
peaks on lattice size. 

To examine whether the incommensurate peak along 
the diagonal direction is produced by a finite 0-0 hop- 
ping tpp, in Fig. 9 (a) and (b) we display the spin struc- 
ture factor S{k) as a function of tpp for different e. The 
parameters are the same as in Fig. 7. Here a nonzero tpp 
makes the hole filling correspond to a closed-shell case. 
From Fig. 9 (b), even for a very small tpp, a weak peak 
at the (tt — S',Tr — 6') still exists. With increasing tpp, 
the spin-spin correlations are strongly suppressed near 
the AF wavevector (tt, tt), and at the same time the am- 
plitude of the incommensurate peak along the diagonal 
direction or tendency to this peak forming is enhanced. 
For the half-filling case (data not shown) we also observed 
that increasing tpp greatly suppresses AF order. 



Finally, we report how the AF long range order in 
the half-filling case is destroyed by the hole doping. In 
Fig. 10(a) and (b), the spin structure factor is plotted for 
different hole filling cases. From Fig. 10(a), it is clearly 
seen that AF spin-spin correlation is strongly suppressed 
by the hole doping. As shown in Fig. 10(b), in the light 
hole doping region {Nh — 42), there exist two incommen- 
surate peaks. When the system is doped to Nh = 46, the 
incommensurate peak along the diagonal direction disap- 
pears, but the peak on the face of of the Brillouin zone is 
robust. In the heavy doping region {N^ = 54), the spin 
structure factor is featureless near (7r,7r), and the peak 
occurs at (tt, 0). 



C. Pairing Correlation Functions 

A typical pairing correlation function as a function of 
distance and Vpd is shown in Fig. 11. Here we show the 
d-wave function and see that it is dominated by a large 
peak at short distances (i? < 2). At these distances, in- 
creasing Vpd increases the magnitudes of the correlations 
slightly. At larger distances {R > 2), the trend reverses. 
The dominance of the local peak is such that a measure 
of pairing like the k = dependence of the Fourier trans- 
form of the pairing correlation function or the integral of 
P{R) with a large distance cut-off can exhibit behavior 
indicative of the only the short-range behavior and hide 
the more relevant long-range behavior. 

The short-range behavior is illustrated in more detail 
in Fig. 12. For two different values of Ud, we plot the 
value of the R — peak in the vertex contribution for 
both extended s-wave and d-wave symmetries for several 
values of e. In both the charge-transfer and mixed valent 
regimes, the R = value for both symmetries increases 
monotonically with increasing Vpd- 

For comparison, the long-range {R > 2) behavior is 
illustrated in Fig. 13. In both the s and d-wave chan- 
nels, the long-range vertex contributions in the charge- 
transfer regime decrease with increasing Vpd but in the 
mixed-valent regime it increases. Over the range of Vpd 
simulated, the long-range part of the d-wave contribution 
is consistently larger than the s-wave contribution. 

As a function of filling, the behavior is more complex. 
In Fig. 14 is the short-range part of the vertex contribu- 
tion. In the mixed valent regime, both the extended s- 
wave and d-wave channels decrease rapidly with electron 
and hole doping. In the charge-transfer regime, the d- 
wave functions falls with doping but the extended s-wave 
contribution initially increases. The long-range contribu- 
tion, shown in Fig 15, shows the dominance of the d-wave 
channel. In the mixed-valent regime it basically decreases 
with doping, whereas in the charge-transfer regime, it ini- 
tially increases but only to decrease rapidly for large hole 
doping. Our findings suggest that the conclusions of Dopf 
et al.El reflect the behavior of local interaction vertex, not 
the long range property. 



V. SUMMARY AND CONCLUSIONS 

We summarize our results as follows: Using the CPMC 
method, we simulated the two-dimensional three-band 
Hubbard model to study its charge, spin, and pairing 
correlations as a function of electron and hole doping, 
and the charge-transfer energy e and the Coulomb re- 
pulsion Vpd between charges on neighboring Cu and O 
lattice sites. We found that increasing Vpd suppressed 
the charge-charge correlations between neighboring Cu 
and O sites. In the mixed-valent regime it had the effect 
moving small amounts of charge from the Cu sites to the 
O sites. In the charge-transfer regime, the effect was the 
opposite. Upon hole doping, more of the extra holes went 
to the O sites than to the Cu sites. 

At a given doping, increasing Vpd increased the spin- 
spin correlations in the charge-transfer regime but de- 
creased them in the mixed-valent regime. Also electron 
and hole doping away from half-filling was accompanied 
by a rapid suppression of anti-ferromagnetic correlations. 
As a function of doping, e, and Vpd, the behavior of the 
magnetic moment on the Cu sites was strongly correlated 
with the behavior of the charge on the Cu sites. 

As a function of distance, both the d^2_y2--wave and 
extended s-wavc pairing correlations decayed quickly. In 
the charge-transfer regime, increasing Vpd decreased the 
long-range part of the correlation functions in both chan- 
nels, while in the mixed-valent regime, it increased the 
long-range part of the s-wave behavior but decreased 
that of the d-wave behavior decreased. Still the d- 
wave behavior dominated. At a given doping, increasing 
Vpd increased the spin-spin correlations in the charge- 
transfer regime but decreased them in the mixed-valent 
regime. Also electron and hole doping away from half- 
filling was accompanied by a rapid suppression of anti- 
ferromagnetic correlations. 

We presented a more extensive study of the effects of 
Vpd than previous QMC studies at zero and finite temper- 
ature. Our results illustrate the presence of both s and 
d-wave correlations in the charge-transfer regime, with 
the d-wave correlations generally dominating. These re- 
sults highlight the difference in the behavior between the 
short and long range part of these correlation function. 
Figures of merits that include the short-range part are 
dominated by the behavior of the short-range part. For 
the system size studied the correlations are weak. We 
did not study these correlations as function of system 
size. We believe the size dependence will be the same as 
previous studies. 

Lastly, as a complementary part to results presented 
so far, we briefly report the effects of Up on charge, mag- 
netic and pairing correlations. In both the mixed valent 
and charge-transfer regions, and for all band fillings, we 
found that a finite Up {Up = 2.0) moves some charge 
from the oxygen sites to the copper sites, which causes 
an increasing of magnetic moment at the copper sites. 
Consistent with previous observationsE3E3 that Up has a 



negative effect on the hole binding, the long-range part of 
the d-wave correlation functions is suppressed by Up. Our 
simulation results also show that Up has a larger effect in 
the mixed valent region than that in the charge-transfer 
region. 
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FIG. 1. Phase convention for the hopping matrix elements. 
The copper d^2_y2 orbital is surrounded by the oxygen px and 
Py orbitals. The hopping matrix elements are shown with 
their corresponding phase. 



FIG. 2. The band-structure of an infinite systems for e 



and tpp/tpd = 0.3. 



FIG. 3. Fermi surfaces of an infinite-system for e = 3 and 
tpp/tpd ~ 0.3. From the inside out, the hole fillings are 54/36. 
46/36, 42/36, 36/36 (dashed line), and 26/36. 



FIG. 4. Average charge on Cu sites as a function of Vpd for 
different fillings and charge-transfer energies e. 



FIG. 5. Charge correlation between neighboring Cu and 
O sites as a function of Vpd for different fillings and 
charge-transfer energies e. 



FIG. 6. Average magnetic moment at the Cu sites as a 
function of Vpd for different fillings and charge-transfer ener- 
gies e. 



FIG. 7. Average static spin structure factor for Cu sites as 
a function of the wavevector k and Vpd- e = 1 and the number 
of holes equals 42. 



FIG. 8. Average static spin structure factor for Cu sites as 
a function of the wavevector k and Vpd- e = 3 and the number 
of holes equals 42. 



FIG. 9. Average static spin structure factor for Cu sites as 
a function of the wavevector k and tpp. (a) e = 1, (b) e = 3, 
and the number of holes equals 42. 



FIG. 10. Average static spin structure factor for Cu sites 
as a function of the wavevector k and A^'^ with e = 3. (a) 
Half-filling and doped cases and (b) Doped cases. 



FIG. 11. d-wave pairing correlation function as a function 
of distance R for different values of Vpd- The number of holes 
equals 42. 



FIG. 12. Local (R = 0) vertex contributions to the ex- 
tended s and d-wave pairing correlation function as a function 
of Vpd for different values of e and Vpd ■ The number of holes 
equals 42. 



FIG. 13. Long-range (averaged for R > 2) part of the ver- 
tex contributions (averaged over i? > 2) to the extended s 
and d-wave pairing correlation function as a function of Vpd 
for different values of e and Vpd- The number of holes equals 
42. 



FIG. 14. Local part {R — 0) of the vertex contributions to 
the extended s and d-wave pairing correlation function as a 
function of filling. Vpd = O.Thc number of holes equals 42. 
(a) e = 1. (b) e = 3. 



FIG. 15. Long-range part (averaged for R > 2) of the ver- 
tex contributions to the extended s and d-wave pairing corre- 
lation function as a function of filling. Vpd ~ O.The number 
of holes equals 42. (a) e = 1. (b) e = 3. 
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